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Versions of GMRES with deflation of eigenvalues are applied to lattice QCD problems. Approximate eigenvec- 
tors corresponding to the smallest eigenvalues are generated at the same time that linear equations are solved. 
The eigenvectors improve convergence for the linear equations, and they help solve other right-hand sides. 



1. Introduction 

This paper looks at the iterative solution of 
complex, non-Hermitian systems of linear equa- 
tions associated with a lattice QCD problem in 
particle physics. Let the n by n system of equa- 
tions be Ax — b. We use a version of the GM- 
RES method that deflates eigenvalues and can 
solve systems with multiple right-hand sides. We 
give examples that show deflating eigenvalues can 
make a significant improvement in the conver- 
gence for QCD matrices. 

We work with the non-Hermitian Wilson-Dirac 
matrix M(x, y) given by, 

M(x, y) = I -«£)[(! + -/M(x - a M ) 



Sx,y+a„ + (1 - lnW^XjS^y-aJ, 



(1) 



because we perfer to consider systems for which 
multiple-mass shifts are possible. Our lattice here 
is of size 12 3 x 24. After the standard even/odd 
preconditioning, the lattice QCD problem has the 
form, 



A e x e — b e + nD eo b , A e — / k D eo D 



(2) 



where A e is of dimension 248,832. We set k = 
0.1575, which for this size matrix and beta value 
(6.0) means working essentially at 

^critical • 

We 

consider a typical, non-exceptional configuration, 
Ufj,. The method is designed to efficiently solve 
the multiple right-hand sides of Ax = b that oc- 
cur when the "all to all" propagators used for dis- 
connected diagrams are calculated. We will use 
multiple right-hand sides formed with Z(2) noise 



vectors. The matrix is not extremely sparse, with 
about 200 nonzeros per row. 

Krylov subspace methods are iterative meth- 
ods for solving large systems of linear equations. 
With approximate solution xq, and residual vec- 
tor ro = b — Axq, the linear equations can be 
recast as A(x — xq) — r^. Then the Krylov 
subspace of dimension m for this problem is 
Span{ro, Aro, A 2 r , A m ~ 1 r }. The conjugate 
gradient method is for Hermitian problems, and 
GMRES g|| is a well-known Krylov method for 
the non-Hermitian case. 

The effectiveness of Krylov methods is of- 
ten controlled by the distribution of eigenval- 
ues. The existence of small eigenvalues can sig- 
nificantly slow the convergence rate. This moti- 
vates attempts to remove or deflate some eigen- 
values from the effective spectrum for an iterative 
method. 



2. The Method 

The GMRES method fully orthogonalizes a ba- 
sis for the Krylov subspace. As the iteration pro- 
ceeds, the expense and storage both grow. This 
makes occasional restarting necessary. We refer 
to each pass through the GMRES iterations be- 
tween restarts as a "cycle" . 

Several approaches have been proposed for de- 
flating eigenvalues from GMRES; see || for refer- 
ences. For deflation in QCD problems, see j| (/J. 
Our approach differs from previous QCD defla- 
tion in that the eigenvalue problem is solved si- 
multaneously with the linear equations. 
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In this paragraph, we describe the method for 
the first right-hand side. After the first GM- 
RES cycle, approximate eigenvectors called har- 
monic Ritz vectors, are computed from the same 
Krylov subspace generated for solving the linear 
equations. Let the approximate eigenvectors cor- 
responding to the k approximate eigenvalues of 
smallest modulus be yi,y%, ■ . . j/fc. For the next 
cycle, the subspace used for GMRES is 

Span{r , Ar , A 2 r ,A 3 r Q , A™-*- 1 ^, 

yi,lk,...yk}- (3) 

Having the approximate eigenvectors in the sub- 
space accomplishes two things. The correspond- 
ing eigenvalues are deflated once the approxi- 
mate eigenvectors are moderately accurate. Also 
the approximate eigenvectors are improved as 
the method proceeds. It sometimes takes sev- 
eral cycles to improve the approximate eigen- 
vectors to the point that they are useful. This 
method is called GMRES-DR (GMRES with de- 
flated restarting) g. GMRES-DR gives the ap- 
proximate eigenvectors in the form of a short 
Arnoldi-type recurrence 

AV k = V k+1 H k , (4) 

where Vk is an n by k orthonormal matrix whos 
columns span the subspace of approximate eigen- 
vectors, and Hk is a full k + 1 by k matrix. 

Now we look at solving the second and sub- 
sequent right-hand sides. We use the eigenvec- 
tor information that was generated in solving the 
first right-hand side. An efficient approach in 
which the eigenvalues are deflated outside of the 
GMRES cycles with a simple projection is called 
GMRES-Proj @. With the short Arnoldi-type 
recurrence (^) from GMRES-DR, we need to store 
only fc + 1 vectors of length n in order to have 
access to both the approximate eigenvectors and 
their products with A. This allows for fairly inex- 
pensive projections. We use a minimum residual 
projection || here. The GMRES-Proj method 
applies cycles of standard GMRES, with a projec- 
tion over the approximate eigenvector subspace 
in between cycles. This projection is not needed 
between all of the GMRES cycles. 
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Figure 1. Solution of First Right-hand Side 



3. Experiments 

We test the approaches just discussed with the 
matrix A e , given earlier. 

Example 1. For the first right-hand side, we 
compare GMRES-DR with standard GMRES. 
Figure |l| gives a plot of the residual norms for 
GMRES-DR(40,20), GMRES-DR(20,10), GM- 
RES(40), and GMRES(20). GMRES-DR(40,20) 
uses subspaces of dimension 40, of which 20 ba- 
sis vectors are approximate eigenvectors. GM- 
RES(40) also restarts when the subspace reaches 
dimension 40. Starting around iteration 300, 
GMRES-DR(40,20) performs considerably better 
than the other methods. At that point, it has de- 
veloped good enough approximations to some of 
the smallest eigenvalues. The eigenvalues halfway 
surround the origin ||, on the positive real side. 
This makes the problem difficult, but removing 
some of the surrounding eigenvalues nearest the 
origin is helpful. 

Example 2. For the second right-hand side, 
we compare GMRES-Proj with standard GMRES 
and also with GMRES-DR. With GMRES, there 
is no deflation, and with GMRES-DR, the defla- 
tion happens only after accurate enough approxi- 
mate eigenvectors develop. Meanwhile, GMRES- 
Proj is able to deflate from the start. Fig- 
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4. Conclusion 

Deflating eigenvalues is useful for lattice QCD 
problems, particularly for the second and subse- 
quent right-hand sides. While Z(2) noise vec- 
tors were used for these examples, the speedup 
from deflation is independent of the nature of the 
right-hand side. Future plans include implement- 
ing multiple mass versions of GMRES-DR and 
GMRES-Proj that solve several systems of equa- 
tions with shifted matrices but the same right- 
hand side. We also would like to investigate de- 
flating eigenvalues for the subsequent right-hand 
sides from Lanczos methods such as the conjugate 
gradient method and BiCGSTAB. 



Figure 2. Solution of Second Right-hand Side 



ure H has the residual norms for GMRES(20)- 
Proj(30), GMRES(20)-Proj(20), GMRES(20)- 
Proj(lO), GMRES(20), and GMRES-DR(40,20). 
GMRES(20)-Proj(30) refers to cycles of GM- 
RES(20) with occational projections over 30 ap- 
proximate eigenvectors in between. We project 
in between every third cycle. GMRES(20)- 
Proj(20) and GMRES(20)-Proj(10) use the ap- 
proximate eigenvectors generated while solving 
the first right-hand side with 400 iterations of 
GMRES-DR(40,20) and GMRES-DR(20,10), re- 
spectively. GMRES(20)-Proj(30)'s eigenvectors 
come from 610 iterations of GMRES-DR(50,30) 
(more iterations are needed in this case for the 
eigenvectors to become accurate enough). The 
expense per iteration is nearly the same for 
GMRES(20) and all the GMRES-Proj methods, 
while GMRES-DR(40,20) has greater orthogo- 
nalization expense. More storage is needed for 
the GMRES-Proj methods than for GMRES. 
For example, GMRES (20)-Proj (20) needs over 
40 vectors of length n, while GMRES(20) uses 
a little over 20. The deflation in GMRES- 
Proj makes it much better than GMRES (20). 
GMRES (20)-Proj (20) is also considerably better 
than GMRES-DR(40,20), because the deflation 
can start from the beginning (both use 20 Krylov 
vectors and 20 eigenvectors for each cycle). 
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